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We present a numerical study of thermodynamical properties of dimerized frustrated Heisenberg 
chains down to extremely low temperatures with applications to CuGe03 . A variant of the finite 
temperature density matrix renormalization group (DMRG) allows the study of the dimerized 
phase previously unaccessible to ab initio calculations. We investigate static dimerized systems as 
well as the instability of the quantum chain towards lattice dimerization. The crossover from a 
quadratic response in the free energy to the distortion field at finite temperature to nonanalytic 
behavior at zero temperature is studied quantitatively. Various physical quantities are derived and 
compared with experimental data for CuGe03 such as magnetic dimerization, critical temperature, 
susceptibility and entropy. 
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I. INTRODUCTION 

The investigation of low dimensional quantum spin sys- 
tems has attracted widespread and general interest. Ac- 
tive research is performed on a large class of magnetic 
materials, experimentally as well as theoretically. Up 
to date, an enormous understanding of the low-energy 
physics of quasi-one dimensional systems such as spin- 
Peierls, Haldane and ladder systems has been reached. 
Still it is a challenging task to investigate especially the 
low temperature regime of these systems, because for T 
close to zero thermal fluctuations are strongly reduced 
and quantum fluctuations become dominant. 

The recently developed transfer-matrix DMRG, which 
combines White's DMRG idea |IJ] with the transfer- 
matrix approach [^) is particularly suited for this task 
[^-^|, because the thermodynamic limit in quantum 
space can be performed exactly. Most other numer- 
ical techniques such as the "traditional" Hamiltonian 
DMRG or exact diagonalization are restricted to either 
zero temperature or to rather large temperatures. Quan- 
tum Monte Carlo algorithms f| also can only be applied 
to systems finite in both quantum and Trotter direction. 
Furthermore, the sign problem for frustrated systems is 
a rather sophisticated point associated with the latter 
approach. 

In this paper we investigate the thermodynamics of 
dimerized frustrated Heisenberg chains, which are be- 
lieved to model appropriately the spin degrees of free- 
dom of, for instance, the inorganic spin-Peierls compound 
CuGe03 . The purpose of this paper is twofold. First, 
we want to show that a density matrix renormalization 
group treatment of a "lattice path integral formulation" 
for general spin chains with dimerization and frustra- 
tion is possible. Within this approach we will under- 
stand qualitatively and quantitatively the instability of 
one-dimensional frustrated antiferromagnetic Heisenberg 
systems towards dimerization at finite temperature. Sec- 
ond, we want to show that many of the thermodynamic 



properties of CuGe03 can be understood within the one- 
dimensional adiabatic approach. At least with respect 
to the location of T SP , the entropy and the susceptibility 
down to temperatures of 5K the role of higher dimen- 
sional magnetic couplings is less relevant. 

The paper is organized as follows. Computational as- 
pects and the model Hamiltonian are given in section 
II. In section III we analyze the dependence of the free 
energy on the dimerization and the crossover from ana- 
lytical to non-analytical behavior in the limit of vanishing 
temperature. We compare our numerical findings to the- 
oretical results based on scaling arguments. In addition 
we present the temperature dependence of the order pa- 
rameter within the adiabatic approach. In section IV we 
show results for the magnetic susceptibility as well as the 
entropy which compare excellently with the experimental 
data for CuGe03 . Summary and conclusion are given in 
section V. 



II. MODEL AND METHOD 

The dominant magnetic interactions in CuGe03 are 
due to Heisenberg spin exchange between Cu 2+ ions 
along the c-axis of the crystal with strong nearest (J) and 
next-nearest neighbor (a J) interactions |t]] . As usual, we 
assume a linear dependence of the nearest neighbor ex- 
change integral on the local structural deformation. The 
commonly used Hamiltonian in adiabatic approximation 
reads 



H = 



J2\ J((1 + <5,)S,S 



i+l 



(1) 



where Si are spin 1/2 operators, <5; = (— l) l <5 denotes the 
modulation of the exchange coupling in the dimerized 
phase (<5 > is assumed) and h = gusH^ represents 
the effective external magnetic field. The last two terms 
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in (|l]) are the elastic energy associated with the lattice 
distortion and the Zeeman energy. 

In order to have a translationally invariant description 
we combine each two spin 1/2 to an effective four state 
object resulting in a model with only nearest neighbor in- 
teraction. A similar strategy has been adopted in Ref. ||] 
to cast a system with next-nearest neighbor interaction 
in a form suitable for applications of the transfer matrix 
DMRG algorithm. Firstly, the elastic part of the Hamil- 
tonian is neglected pending further notice, i.e. K = 0, 
since it only enters as a trivial constant in all energies 
and can be added finally. 

After chequerboard decomposition B of the effective 
four state problem the partition function is calculated 
by means of the quantum transfer matrix Tm P,|l0|. 
In the thermodynamic limit all physical quantities are 
determined by the largest eigenvalue A of Tm and the 
corresponding left and right eigenstate. More precisely, 
the free energy per site of the spin- 1/2 chain reads 
/ = — -j In A. Expectation values of local operators like 
the internal energy u, magnetization m or nearest neigh- 
bor correlation (S^Si+i) etc. can be expressed in terms 
of the eigenstates (see Q and for a general reference 
to transfer matrix formalisms e.g. ]i"l"|]). The quantum 
transfer matrix Tm, A and the eigenvectors are calculated 
by iterative augmentation in Trotter direction within an 
application of the DMRG concept as in f§§. For all cal- 
culations we kept 24 states in the renormalization step 
and chose (TM) -1 = 0.05. 



The relation a(T) = K - JA(T) holds exactly, but we 
are restricted to finite 8 in our numerics since the error 
in A(T, 8) scales with 1/8. (Determining A(T) directly 
from /(T, 8) — f(T, 0) would even involve an error of or- 
der 1/8 2 .) In the sense of a multi-point formula we can 
assume an expansion like (Q) up to order 8 2n to hold ex- 
actly, and calculate the corresponding n coefficients by 
considering the function A(T, 8) for n different values of 
8. In this way we obtain reliable approximations for the 
functions a(T), b(T), .... In the present work we employed 
a n = 3 point formula and used the three smallest dimer- 
izations accessible. 
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III. GENERAL RESULTS 

First we calculate the free energy /(T, 8) for the system 
(|l|) for a series of fixed distortions Si = (— From 
general principles we expect analytic behavior for any 
finite temperature, i.e. 



f(T,8)=f(T,0) 



l a (T)8 2 + ±b(T)8* 



0(8% (2) 



where the odd powers in 8 disappear due to the obvi- 
ous symmetry (8 — > —5). The coefficient a(T) of S 2 is 
the reciprocal of the static "dimerization susceptibility" . 
Hence, a zero in this coefficient signals a structural tran- 
sition. 

Numerically we determine the coefficient by analyzing 



A(T,6) = j 



l df(T,6) 
8 88 



K 



(giga) - (S 2 S 3 ) 
28 



(3) 



which in the limit 8 — > is simply denoted by A(T). Note 
that A is independent of K. In Fig.l we show DMRG 
results for the temperature dependence of A(T, 8) for var- 
ious (5 in the case a = 0. 



FIG. 1. A(T, 8) as defined in the text versus temperature 
for a = and 8 = 0.15, 0.1, 0.05, 0.04, 0.03, 0.01 (from bottom 
to top) and the extrapolated A(T) (dashed line). The inset 
shows a log-log plot of (A(T) + c) with c = —1.14 (dashed 
line) together with a line with slope -1 (thin solid line). 

Here we are particularly interested in A(T) at low tem- 
perature showing a divergence at T = 0. This singular- 
ity is caused by the failure of (^|) at zero temperature. 
From jl2],[l3| and also from scaling relations of the quan- 
tum critical Heisenberg chain flifl we know for a < a c 
(a c w 0.2412 @,0) that f(0,8) - /(0,0) ~ 8 4/3 . For 
a > a c the quantum chain shows spontaneous magnetic 
dimerization as well as an excitation gap which leads to 
/(0,S)-/(0,0)~<5. 

More quantitatively, we first define a temperature vari- 
able t = f(T,0) - /(0,0). The singular part of the 
free energy in dependence of t and 8 is denoted by 
f(t, 8) = f(T, 8) - f(0, 0). Then, we expect the following 
scaling relation 



f(t,8) 



)f(ti x ,siy). 



(4) 



(See e.g. |Tj| as a general reference.) The exponents x 
and y are to be determined. From (0) we derive in the 
usual manner f{t,8) = t^ x f(l, S/t^). Fixing 8 = 
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directly implies x — 1 due to the definition of t. Next we 
use the fact that f(Q,S) ~ 6 V which implies f{l,z) — z v 
for large z and furthermore y — I/77. Lastly we note 



d 2 - Id 2 



/(l,0) -('-;. (5) 



Fig. [I] displays A(T, 5) for a = and various 8 and 
A(T) (dashed line) derived from the curves with the low- 
est dimerization. As expected the low-temperature be- 
havior of the free energy for vanishing dimerization is 
found to be t = T 2 with an extrapolated groundstate en- 
ergy fo — —0.443. Here we have rj = 4/3, thus we should 
observe A(T) ~ l/T at sufficiently low T. This is ex- 
cellently confirmed by our numerics (see inset of Fig. [I]). 
As the scaling relations only apply to the divergent terms 
we are allowed to subtract regular terms. For instance an 
additive constant for A(T) does not change the scaling, 
but enlarges the temperature range where this sets in. 

For overcritical frustration a > a c a similar scaling 
analysis can be performed. The main difference is the 
absence of quantum criticality due to the gap A of the 
elementary excitations. 
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FIG. 2. A(T, S) versus temperature for a — 0.5 and 
5 = 0.02, 0.015, 0.01, 0.005, 0.003 (from bottom to top) as well 
as the extrapolated A(T) (dashed line). The inset shows 
ln(A(T,S)(T/J) 3/2 ) + cT/J with c = 0.56 (solid lines) ver- 
sus J/T, the same transformation of the extrapolated A(T) 
(dashed line) and a line with slope A ST /2 (thin solid line). 

Fig. ^ shows A(T, 8) for a = 0.5 and various 8. 
The coefficient A(T) (dashed line) is deduced as before. 
The temperature variable is t = T 3 / 2 exp(— A/T) where 
A = A ST /2 is half the singlet-triplet gap for 5 = 0. 
The origin of this factor 1/2 can be understood. The 
elementary excitations are spinons each with a gap A, 
however occurring pairwise. Therefore the spectroscopic 
data show a gap 2 A = A ST . For the usual thermody- 
namical data like the free energy the activated behavior 



shows a characteristic energy corresponding to the gap 
of the individual elementary excitations whether these 
occur in arbitrary numbers or restricted to even num- 
bers does not matter, see e.g. |Q. (A different situ- 
ation is realized for 8 > due to the confinement of 
spinons in pairs. Here a crossover in the activated be- 
havior sets in at temperatures of the order of the binding 
energy E B «-» T B . At temperatures above T B the gap 
in the activated behavior /2, below T B it is A ST .) 

For 8 = the fitted singlet-triplet gap A ST » 0.23(2) 
compares fairly well to the value 0.2338(6) jl9|,^0| for 
the Majumdar- Ghosh model. The extrapolated ground- 
state energy reads fo = —0.3749, also close to the ex- 
act result — |. Here we have to set 77 = 1, which im- 
plies A(T) ~ t . Indeed we find a linear regime for 
\n[A{T){T / Jf/ 2 ] as a function of l/T (see inset of Fig. 

t, which confirms the prediction of the scaling ansatz 
). Again, the subtraction of regular terms is admis- 
sible. The slight deviations from slope A ST /2 must be 
apportioned to an inaccuracy in the extrapolated A(T). 

Up to now, we have considered static dimerized sys- 
tems. If (5 is treated as a thermodynamical degree of 
freedom, the system favors the dimerization which min- 
imizes the free energy. The spontaneous dimerization 
8(T) as function of temperature can either be calculated 
by comparing the free energies for different dimerizations 
S (K > fixed) or due to (||), S(T) is implicitly defined 
by 



A(T,5(T)) = j, 



(6) 



i.e. the intersection of A(T, 6) for fixed 5 with the con- 
stant K/J. The transition temperature T S p is defined 
by <5(r s ~) = 0. For higher (lower) temperatures the free 
energy is increased (lowered) for non-vanishing S corre- 
sponding to the uniform (dimerized) phase. For deter- 
mining the spontaneous dimerization we use (||) rather 
than the free energy intersection whose accuracy is ex- 
pected to be a factor 1/5 worse. 

In Fig. H the temperature dependence of the spon- 
taneous dimerization (S(T)) is depicted exemplarily for 
a = 0, K = 1.6 J and a = 0.35, K = 2.4J. The exponents 
at the critical temperature T SP are close to 1/2 which is 
expected due to the mean-field treatment of the elastic 
degrees of freedom. The numerical analysis has been per- 
formed for a discrete set of parameters 8 with maximum 
value 0.15 which appears to be an upper bound for the 
actual 8{T = 0) of both examples. The value is in agree- 
ment with earlier self-consistent calculations by means of 
the DMRG for zero temperature pl[ . 

Moreover, the critical temperature as a function of the 
elastic constant as well as a function of the dimerization 
S(T = 0) can be deduced from A(T, 8). The transition 
temperature of CuGeO-3 will be discussed in the following 
section. 
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FIG. 3. Spontaneous dimerization versus temperature for 
a = 0, K = 1.6J (circles) and a = 0.35, K = 2.4J (dia- 
monds). The lines show fitted functions of type c ■ (T c — T) 
with (3 fixed at 1/2. 

IV. APPLICATION TO CuGe0 3 

Since the discovery of CuGeC>3 as the first inorganic 
spin-Peierls compound there have been numerous stud- 
ies of appropriate sets of microscopic couplings |l7],|22|,|3| . 
Within the one dimensional adiabatic approach com- 
monly used parameter are: J « 160-ftT, a ks 0.35 and 
5(T = 0) « 1.3% These parameters were 

derived mainly by fitting the experimental susceptibil- 
ity data in the uniform phase well above T SP and from 
the singlet-triplet gap (singlet-singlet gap p4| ) at zero 
temperature. Note, that the entropy for these interac- 
tion parameters also satisfies bounds derived from the 
specific heat data p3| |. 

Next we focus on a brief comparison to experimental 
data. Fig. [I] shows experimental and theoretical suscep- 
tibility data. The agreement for T > 35K has already 
been proven in p3[ . We also find agreement in the re- 
gion between T SP and 35K. The deviation close above 
the phase transition may be interpreted as a precursor 
due to fluctuations in the real system while our theoret- 
ical calculations indicate a sharp transition point. Using 
an elastic constant of A' = 11J the calculated critical 
temperature is equal to the experimentally determined 
value T S p = 14. AK. Below T SP the susceptibility is sys- 
tematically overestimated. With K = 10.2 J, the DMRG 
results perfectly coincide with experiment (see Fig. 
but the transition takes place at about 15.2K. Of course, 
this a natural effect due to the meanfield approximation. 
A similar overestimation of T SP is observed by adjusting 
a square root behavior to the experimentally measured 
lattice dimerization |p5f . 

In Fig. | the entropy S(T,5(T)) for a = 0.35 is de- 



picted. The consistency with the entropy bounds derived 
by an analysis of the specific heat was already mentioned 
in p!| for temperatures above 35K. The experimentally 
determined bounds are also respected for temperatures 
between T SP and 35-RT, below the critical temperature we 
observe deviations of about 5% with respect to the up- 
per bound for K = 11J. Decreasing the elastic constant 
to K — 10. 2J with the consequences stated above the 
bounds are respected in the entire temperature range. 
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FIG. 4. Zero field susceptibility of CuGe03 . Comparison 
of experimental data (solid line) for the c-axis (p-factor 2.05) 
with exact diagonalization (L = 16) of the non-dimerized frus- 
trated Heisenberg chain (dashed line) and DMRG results us- 
ing K — UJ (circles), K — 10. 2J (crosses). 

Now the low temperature regime is accessible for the 
first time by means of the transfer matrix DMRG. To 
derive the experimentally observed thermodynamics we 
have to set K = 11J (10. 2 J) which implies a maximum 
dimerization of 5(T = 0) w 2.6% (2.8%) and A ST w AOK 
(42 K). The larger dimerization compares well with re- 
cent results of Biichner et al. [ p6[ | . They conclude a min- 
imum dimerization of 3% based on measurements of the 
pressure dependence of the exchange couplings together 
with the structural distortion. Interestingly, the larger 
saturation dimerization leads to a singlet-triplet gap close 
to the average gap Ast ~ 44 K 1 which may be used within 
the one dimensional approach regardingthe interchain 
interaction (b direction) approximately ]2jfl . 

Nevertheless, these elastic constants deviate from ear- 
lier T = calculations, where K = 18J (S(T = 0) « 
1.4%) is found to describe the zero temperature physics 
correctly [ p8| . In particular, the experimentally observed 
gap and the critical magnetic field at the D/I transition 
are reproduced for K = 18J. We have to use a smaller 
value for K to obtain accordance with the experimental 
transition temperature. The apparent discrepancy of our 
finite T calculations and the results in (28| may be ex- 
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plained by residual interchain interactions as discussed 
above. 




FIG. 5. Entropy for CuGeOs . Comparison of experimen- 
tally determined bounds (solid lines) with DMRG results us- 
ing K = 11J (circles) and K = 10. 2 J (crosses). 



V. SUMMARY 

By applying a quantum transfer matrix DMRG al- 
gorithm to the frustrated Heisenberg chain with fixed 
dimerization we calculated physical quantities down to 
temperatures of less than 0.05 in units of the coupling 
constant J. 

The results compare to analytical predictions based on 
scaling arguments with perfect agreement. In particular, 
we can clearly distinguish between the quantum criti- 
cal case (a < a c ) and overcritical frustration (a > a c ) 
which are characterized by different divergence behav- 
ior of the ^-coefficient in the free energy. Furthermore, 
the instability towards dimerization can be understood 
within this numerical approach. It is the first calculation 
of the spin-Peierls temperature in the adiabatic model 
(|I|) for realistic sets of parameters J, a, K and the first 
quantitative investigation of the dimcrized phase at finite 
temperatures. 

Especially with respect to CuGeOs , taking the param- 
eters J = WOK, a = 0.35 established by high tempera- 
ture studies and an elastic constant K = 11J (10. 2J), a 
critical temperature of T s t p eory « 14AK (\h.2K) is found. 
Results from the latter parameter set are in perfect agree- 
ment with the experimental results, e.g. magnetic sus- 
ceptibility (Fig. ||) and entropy (Fig. ||). The slight de- 
viation in the value of T SP can be probably traced back 
to the static treatment of the lattice dimerization which 
overestimates transition temperatures. 

In general, our results extend the numerical analysis to 
much lower temperatures, including the dimerized phase. 
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